################################################################################################
rm(list=ls()) 
x=c("lfe","tidyr","broom","ggplot2")

lapply(x, library, character.only=TRUE) # load the required packages
library(dplyr)
library(starpolishr)
library(data.table)


setwd("Derived Data")


dfdist=fread("distanalysisleaid.csv")

dfdiff=fread("dfbartik2005_40km.csv")
tempdf=read.csv("weatherdf.csv",header=TRUE,stringsAsFactors = FALSE)


state_div=read.csv("../Raw Data/4. ACS Data/us census bureau regions and divisions.csv", stringsAsFactors = FALSE)
state_div$stateabb=state_div$State.Code
dfdiff=merge(dfdiff,state_div,by="stateabb")
dfdiff=merge(dfdiff,tempdf,by=c("leaid","year"),all.x=TRUE)
dfdiff$totalprod40=dfdiff$ayproduction2005.coal_40km+dfdiff$ayproduction2005.gas_40km+dfdiff$ayproduction2005.oil_40km+dfdiff$ayproduction2005.renewable_40km

dfdiff$sharecoal=dfdiff$ayproduction2005.coal_40km/dfdiff$totalprod40
dfdiff$sharegas=dfdiff$ayproduction2005.gas_40km/dfdiff$totalprod40


med2005=mean(dfdiff$sharecoal,na.rm=TRUE)
med2005gas=mean(dfdiff$sharegas,na.rm=TRUE)
med2005=quantile(dfdiff$sharecoal,0.5,na.rm = TRUE)
med2005gas=quantile(dfdiff$sharegas,0.5,na.rm=TRUE)
dfdiff$highind=ifelse(dfdiff$sharecoal>med2005, 1,0) # checked
dfdiff$highindgas=ifelse(dfdiff$sharegas>med2005gas, 1,0) # checked

# post 2012 indicator:
dfdiff$post2012=ifelse(dfdiff$year>=2012,1,0);
dfdiff$highpost2005=dfdiff$highind2005*dfdiff$post2012


dfdiff$yearind=factor(dfdiff$year)
dfdiff$yearind <- relevel(dfdiff$yearind, ref = "2011")
tempvarsavg="ay_mean_snow  +aymean_high_num_snow+aymean_severe_num_snow+ay_avg_mean_max + ay_avg_mean_min + ay_avg_bint80  + ay_avg_bint90  + ay_avg_bint100 + ay_avg_bint20 + ay_avg_bint10+ ay_avg_bint0+ ay_mean_prcp + aymean_high_num_prcp+aymean_severe_num_prcp"
demogvar="perfrl + perblk + perhsp + perwht  + malpct_test + pcttested+ totenrl +  perell + perspeced "
econvar="pctemploy + pctmanufac + single_momall + baplusall + pctutility + pctlaborforce"





acsdf=read.csv("acs2009_2018.csv",stringsAsFactors = FALSE, header=TRUE)
acsdf$pctemploy=as.numeric(acsdf$pctemploy); acsdf$pctlaborforce=as.numeric(acsdf$pctlaborforce); acsdf$pctutility=as.numeric(acsdf$pctutility); acsdf$pctmanufac=as.numeric(acsdf$pctmanufac)
acsdf$leaid=acsdf$LEAID; acsdf$LEAID=NULL
dfdiff=merge(dfdiff,acsdf,by=c("leaid","year"),all.x=TRUE)




model1=as.formula(paste("cs_mn_all~ factor(yearind)*factor(highind) + ", demogvar," + ", econvar, "+",tempvarsavg,"  |factor(grade)+ factor(leaid) + factor(subject) |0|leaid",sep="" ))

dfdiffplot=dfdiff[dfdiff$testmonth=="May"  & dfdiff$ayproduction2005.coal_40km>0,]

result1=felm(model1,dfdiffplot,na.action=na.exclude)

tidy_test=tidy(result1,conf.int=TRUE)
coeffkeep=c("factor(yearind)2009:factor(highind)1","factor(yearind)2010:factor(highind)1","factor(yearind)2013:factor(highind)1","factor(yearind)2012:factor(highind)1","factor(yearind)2014:factor(highind)1","factor(yearind)2015:factor(highind)1","factor(yearind)2016:factor(highind)1")

tidy_test=tidy_test %>% filter (term %in% coeffkeep)
yearvec1=c(2009,2010,2012,2013,2014,2015,2016)


tidy_test=cbind(tidy_test,yearvec1)

ggplot(tidy_test,aes(x=yearvec1, y=estimate )) + geom_point() + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high)) + geom_hline(yintercept = 0) + geom_vline(xintercept=2011.5) + geom_point(aes(x=2011, y=0)) + scale_y_continuous(breaks=seq(-0.06,0.06,by=0.02),limits=c(-0.06,0.06)) + scale_x_continuous( breaks = seq(2009, 2016, by = 1),labels=c("08-09","09-10","10-11","11-12","12-13","13-14","14-15","15-16"))+ labs(x="Academic Year", y="Mean Test Score High-Coal vs. Low-Coal Districts") + theme_classic() + theme(axis.text=element_text(size=22),
                                                                                                                                                                                                                                                                                                                                                                                                                                                                     axis.title=element_text(size=22,face="bold"),legend.position="bottom",legend.text = element_text(size = 20),text = element_text(size = 24))
ggsave("../output/Figdiffindiff_coal.pdf")

  